Activity: Spatially Continuous Data I

Activity: Spatially Continuous Data I

Practice questions

Answer the following questions:

1. What is the difference between spatially continuous data and a spatial point pattern?
2. What is the purpose of spatial interpolation?
3. In your own words describe the method of Inverse Distance Weighting.
4. Consider the following spatial interpolation algorithms: Voronoi polygons and k-point means. How do they differ when the number of points used to calculate means is 1?

Learning objectives

In this activity, you will:

1. Explore a dataset with area data using visualization as appropriate.
2. Discuss a process that might explain any pattern observed from the data.
3. Conduct a modeling exercise using appropriate techniques. Justify your modeling decisions.

Suggested reading

• Bailey TC and Gatrell AC [-@Bailey1995] Interactive Spatial Data Analysis, Chapters 5 and 6. Longman: Essex.
• Bivand RS, Pebesma E, and Gomez-Rubio V [-@Bivand2008] Applied Spatial Data Analysis with R, Chapter 8. Springer: New York.
• Brunsdon C and Comber L [-@Brunsdon2015R] An Introduction to R for Spatial Analysis and Mapping, Chapter 6, Sections 6.7 and 6.8. Sage: Los Angeles.
• Isaaks EH and Srivastava RM [-@Isaaks1989applied] An Introduction to Applied Geostatistics, CHAPTERS. Oxford University Press: Oxford.
• O’Sullivan D and Unwin D [-@Osullivan2010] Geographic Information Analysis, 2nd Edition, Chapters 9 and 10. John Wiley & Sons: New Jersey.

Preliminaries

For this activity you will need the following:

• An R markdown notebook version of this document (the source file).

• A package called geog4ga3.

It is good practice to clear the working space to make sure that you do not have extraneous items there when you begin your work. The command in R to clear the workspace is rm (for “remove”), followed by a list of items to be removed. To clear the workspace from all objects, do the following:

rr rm(list = ls())

Note that ls() lists all objects currently on the worspace.

Load the libraries you will use in this activity (load other packages as appropriate).

library(tidyverse)

-- Attaching packages --------------------------------------- tidyverse 1.2.1 --
v ggplot2 3.1.0       v purrr   0.2.5  
v tibble  2.0.1       v dplyr   0.8.0.1
v tidyr   0.8.2       v stringr 1.3.1  
v readr   1.3.1       v forcats 0.3.0  
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(spatstat)

Loading required package: spatstat.data
Loading required package: nlme

Attaching package: 㤼㸱nlme㤼㸲

The following object is masked from 㤼㸱package:dplyr㤼㸲:

    collapse

Loading required package: rpart

spatstat 1.58-2       (nickname: 㤼㸱Not Even Wrong㤼㸲) 
For an introduction to spatstat, type 㤼㸱beginner㤼㸲 


Note: spatstat version 1.58-2 is out of date by more than 3 months; we recommend upgrading to the latest version.

library(spdep)

Loading required package: sp
Loading required package: Matrix

Attaching package: 㤼㸱Matrix㤼㸲

The following object is masked from 㤼㸱package:tidyr㤼㸲:

    expand

Loading required package: spData
To access larger datasets in this package, install the spDataLarge package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`

library(geog4ga3)

Load the data that you will use in this activity:

data("aquifer")

The data is a set of piezometric head (watertable pressure) observations of the Wolfcamp Aquifer in Texas (https://en.wikipedia.org/wiki/Hydraulic_head). Measures of pressure can be used to infer the flow of underground water, since water flows from high to low pressure areas.

These data were collected to evaluate potential flow of contamination related to a high level toxic waste repository in Texas. The Deaf Smith county site in Texas was one of three potential sites proposed for this repository. Beneath the site is a deep brine aquifer known as the Wolfcamp aquifer that may serve as a conduit of contamination leaking from the repository.

The data set consists of 85 georeferenced measurements of piezometric head. Possible applications of interpolation are to determine sites at risk and to quantify uncertainty of the interpolated surface, to evaluate the best locations for monitoring stations.

aquifer$ID <- factor(c(1:nrow(aquifer)))
aquiferID <- as.data.frame(aquifer$ID)

Activity

1. Map the Wolfcamp Aquifer data.

plot_wolfpack <- ggplot(data = aquifer, aes (x = X, y = Y, color = H)) +
  geom_point(alpha = 0.5) + 
  scale_color_distiller(palette = "OrRd", trans = "reverse") +
  coord_equal()
plot_wolfpack

plot_ly(data = aquifer, x = ~X, y = ~Y, z = ~H,
         marker = list(color = ~H, colorscale = c("Orange", "Red"), 
                       showscale = TRUE)) %>% 
  add_markers()

2. Create a surface using Voronoi polygons.

3. Create a surface using IDW. What is the effect of changing the power of the inverse distance function?

W <- owin(xrange = c(0, 200), yrange = c(0, 200))
aquifer.ppp <- as.ppp(X = aquifer, W = W)

39 points were rejected as lying outside the specified window

z_p.idw10 <- idw(aquifer.ppp, power = 10)
z_p.idw2 <- idw(aquifer.ppp, power = 2)
z_p.idw5 <- idw(aquifer.ppp, power = 5)
plot(z_p.idw10)

plot(z_p.idw5)

plot(z_p.idw2)

4. Create a surface using k-point means. What is the effect of changing the number of points used in this algorithm?

target_xy = expand.grid(x = seq(0.5, 100, 2.2), y = seq(0.5, 100, 2.2))
source_xy = cbind(x = aquifer$X, y = aquifer$Y)
kpoint.10 <- kpointmean(source_xy = source_xy, z = aquifer$H, target_xy = target_xy, k =10)
kpoint.1 <- kpointmean(source_xy = source_xy, z = aquifer$H, target_xy = target_xy, k =1)
kpoint.5 <- kpointmean(source_xy = source_xy, z = aquifer$H, target_xy = target_xy, k =5)
ggplot(data = kpoint.10, aes(x = x, y = y, fill = z)) +
  geom_tile() +
  scale_fill_distiller(palette = "OrRd", trans = "reverse") +
  coord_equal()

ggplot(data = kpoint.5, aes(x = x, y = y, fill = z)) +
  geom_tile() +
  scale_fill_distiller(palette = "OrRd", trans = "reverse") +
  coord_equal()

ggplot(data = kpoint.1, aes(x = x, y = y, fill = z)) +
  geom_tile() +
  scale_fill_distiller(palette = "OrRd", trans = "reverse") +
  coord_equal()

5. Discuss the limitations of these approaches. How can you calculate the uncertainty in the predictions?